Coarse-graining dynamics for convection-diffusion of colloids: Taylor dispersion 
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By applying a hybrid Molecular dynamics and mesoscopic simulation technique, we study the 
classic convection-diffusion problem of Taylor dispersion for colloidal discs in confined flow. We 
carefully consider the time and length-scales of the underlying colloidal system. These are, by 
computational necessity, altered in the coarse-grained simulation method, but as long as this is 
carefully managed, the underlying physics can be correctly interpreted. We find that when the 
disc diameter becomes non-negligible compared to the diameter of the pipe, there are important 
corrections to the original Taylor picture. For example, the colloids can flow more rapidly than 
the underlying fluid, and their Taylor dispersion coefficient is decreased. The long-time tails in the 
velocity autocorrelation functions are altered by the Poiseuille flow. Some of the conclusions about 
coarse-graining the dynamics of colloidal suspensions are relevant for a wider range of complex fluids. 

PACS numbers: 



I. INTRODUCTION 

Progress in computer simulations of materials arises 
from two broad classes of innovation. The first is better 
algorithms, for example more efficient ways to calculate 
equations of motion (for molecular dynamics) or to sam- 
ple phase-space (for Monte Carlo). The second comes 
from better coarse-grained models of the underlying ma- 
terials, that is descriptions that are simpler and more 
tractable, but nevertheless retain the fundamental un- 
derlying physics that one is interested in investigating. 
Of course progress also arises simply because computers 
are continually getting faster, so even with old models 
and old methods, new questions can become accessible. 

The focus of this paper will be on new coarse-graining 
methods, especially for non-equilibrium properties of col- 
loidal suspensions. That coarse-graining is needed be- 
comes immediately obvious when considering the enor- 
mous time and length-scale differences between meso- 
scopic colloidal and microscopic solvent molecules. Even 
a nano-colloid of only 10 nm in radius will displace on the 
order of 140,000 water molecules. Moreover, such a small 
colloid would still need about 5 x 10~ 6 s to diffuse over its 
radius, while a typical collision time of water molecules 
is on the order of 10 _15 s. Simulating just a few nano- 
colloids in solution is therefore prohibitively expensive 
with any explicit model of water. On the other hand, it 
is clear from experiments that many general properties of 
colloidal suspensions, such as phase behaviour and basic 
non-equilibrium behaviour, do not depend strongly on 
the radius R of the particles, even though the number of 
water molecules per colloid scales as R 3 . This physical 
fact immediately suggests that coarse-graining methods 
which ignore the fine detail of the solvent should be ap- 
plicable. 



A very fruitful coarse-graining strategy has been 



to 



employ effective potentials that integrate out the solvent, 
and also other degrees of freedom, to create an intuitively 
appealing picture of colloids as "giant atoms" [jj, l2| . Be- 
sides rationalising many results for equilibrium colloids, 
this viewpoint has stimulated the use of colloidal sys- 
tems as a "playground" for probing some of the funda- 
mentals of condensed matter physics. Colloids have the 
advantage that they are much larger than atoms, so that 
they are much more easily visualised. At the same time 
their length-scales are much longer, so that processes in 
the time-domain can be more easily investigated. Ex- 
amples of what can be studied include the kinetics of 
crystallisation Qjthe dynamics of supercooled liquids 0], 
glass formation [1, @ interfacial phenomena [3], disloca- 
tion motion Q and surface melting Q. 

General processes in atomic and molecular materials 
where the non-equilibrium behaviour is dominated by the 
crossing of potential barriers, such as nucleation, glass 
formation and dislocation motion in crystals, are similar 
to the concomitant processes in colloidal systems. But 
at lower colloidal densities, it becomes less obvious that 
this should be the case. This is mainly because, in con- 
trast to atoms, the colloids are in a solvent which induces 
Brownian motion and long-ranged hydrodynamic inter- 
actions (HI) [l(| [H| ■ These HI may decay as slowly as 1 jr 
and can qualitatively affect the dynamical behaviour of 
a suspension. The question of how to coarse-grain away 
the solvent in such a way as to retain both the Brow- 
nian fluctuations and the long-ranged hydrodynamic in- 
teractions is a difficult one. In a previous paper fl2| 
(henceforth paper I) we presented a careful analysis of 
how to achieve this using the stochastic rotation dynam- 
ics (SRD) algorithm of Malevanets and Kapralpjl [T3 |. 
Note that in the literature, this method is sometimes 
also called multiparticlc collision dynamics, see the re- 
cent review by Kapral[15|. SRD has been applied to a 
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wide number of different systems, including fluid vesicles 
in shear flow [161, clay-like colloids [17|, sedimentation of 
colloids (IH 1 1 Oil colloidal rods in shear flow(20l]. knots in 
viral DNA[2ll| and other examples [l5j . 

In the first section of this paper, we briefly review the 
method, and show how it can be used to bridge the time 
and length scales of a colloidal system where HI are im- 
portant. A number of the lessons learned should be valid 
for other coarse-grained studies of dynamics. For exam- 
ple, we argue that many methods effectively use a "tele- 
scoping down" procedure to compress the physical time- 
scales to a more computationally manageable hierarchy. 
Once the simulations are completed, this hierarchy must 
be "telescoped out" in order to make contact with phys- 
ical reality. Of course this cannot be done without losing 
some information; there is no such thing as a free lunch. 

To further illustrate how such a coarse-graining 
method for dynamics can be used to investigate the com- 
bined effects of Brownian motion and HI on colloids, we 
study the diffusive behaviour of colloidal discs confined 
to a narrow two dimensional channels with a fluid under- 
going Poiseuille flow. This convection-diffusion problem 
was first studied in a classic paper by G.I. Taylor [22j| , He 
pointed out that if a solute diffuses between the stream- 
lines of a flow field, this adds an axial diffusive component 
relative to the average flow field, an effect now called 
Taylor dispersion. In his original work, Taylor derived 
the asymptotic form of the dispersion coefficient. Van 
den Broeck [23[ later derived an exact expression for the 
dispersion valid at all times. Rather surprisingly Tay- 
lor found that the effective spreading constant, or Taylor 
diffusion coefficient as we shall refer to it, was inversely 
proportional to the molecular diffusion coefficient D and 
has the form 



where L denotes the channel radius and u the average 
flow velocity for a three-dimensional pipe. Further work 
in the field has lead to the development of a more gen- 
eralised Taylor dispersion theory, as illustrated in the 
review by Van den Broeck [241 ]. Dispersion theory can 
have applications in a multitude of disciplines. Oscil- 
latory flows, such as those that arise in estuaries and 
blood flow have been the subject of intensive study [24J. 
The effect of interparticle interactions, such as chemical 
reactions and different boundary conditions are also of 
interest [24[ . Much of this work has focussed on particles 
whose radius R is much smaller than the channel radius 
a. Here we study in particular what happens to Taylor 
dispersion and related problems when the colloidal ra- 
dius is no longer negligible compared to the pipe radius, 
building on earlier work by Brenner and Edwards [25|. In 
such situations, we find that the colloid flows faster than 
the average of the underlying solvent, mainly because 
the colloid can only sample the faster flow lines near the 
centre of the flow profile. At the same time, the Taylor 
dispersion coefficient is decreased, because the colloid can 



sample a smaller range of flow velocities. For very nar- 
row pipes more subtle hydrodynamic effects due to direct 
interactions with the walls start to become important. 

The paper is organised as follows: In section II we 
discuss some salient points of our coarse-graining scheme 
for colloids. In section III we describe Taylor dispersion 
for point particles, whereas section IV considers the case 
where the colloid radius is no longer negligible compared 
to the pipe width. Finally, in section IV, wc summarise 
our main conclusions. 



II. COARSE-GRAINING COLLOIDAL 
DISPERSIONS 

A. Effective potentials for colloids 

The detailed interactions between colloidal particles 
can be very complex and affected by other (smaller) col- 
loidal particles, polymers, ions and the properties of the 
underlying solvent. However, much progress can be made 
by integrating out these degrees of freedom to generate 
an effective picture with radial potentials [l|, 0]. The ad- 
vantage of these are that the well-developed statistical 
mechanical machinery developed for atomic and molecu- 
lar potentials can now be applied to colloidal suspensions. 
Of course colloids can also be more complicated, for ex- 
ample much recent progress has been made on "colloidal 
molecules" 12611 . and colloids can easily have anisotropic 
interactions 27]. But even for these systems, the picture 
of effective potentials is a powerful abstraction [28j. 

Nevertheless, it must always be kept in mind that these 
effective potentials are not potential energy functions in 
the Hamiltonian sense, but rather include within their 
definition statistical averages over configurations. Their 
character is therefore that of a free-energy. This has a 
number of implications such as transferability problems, 
where a potential derived in one context does not perform 
well in another context. A good example would be if the 
derived potential depends on the overall colloid density. 
Then using it at a different density could lead to errors. 

Another, more subtle, but perhaps equally important 
issue is that of representability problems^^, where the 
effective potential may be used to represent a certain 
subset of physical quantities, but not a different subset. 

For example, consider a one-component colloidal fluid 
interacting with a three-body Hamiltonian of the form: 

H = K + Y j w (2 \n j )+ w (3) (iy,i>,r fcl ), (2) 

i<j i<j<k 

where r\ denotes the position of particle i and = Vi—Yj 
and r>ij — |rj — Tj\. K is the kinetic energy opera- 
tor, w^ 2 \r) is an isotropic pairwise additive potential, 
and «/ 3 ) (rjj , r^, r^) is a triplet or three-body potential. 
Three-body potentials are expensive and cumbersome to 
simulate, and so one might want to coarse-grain them to 
a simpler isotropic pair representation. 
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Henderson [30( | first showed, using arguments very sim- 
ilar to those used by Hohenberg and Kohn [3l[ in their 
famous proof relating the one-body potential to the one- 
body density (which laid the foundation for density func- 
tional theory), that "the pair potential v(r) which gives 
rise to a radial distribution function g(r) is unique up to 
a constant.". An extended proof for orientational corre- 
lations can be found in a book by Gray and Gubbins [12] , 
while a more rigourous mathematical discussion is pro- 
vided by Chayes, Chayes and Lieb (33[. This means that, 
for a given state-point, the g(r) generated by a Hamilto- 
nian like that of Eq. can be reproduced by a unique 
effective pair potential v g (r). We call approaches that 
derive v g (r) the structural route to deriving an effective 
potential. The difference with the bare-pair potential 
u>( 2 ) (r) can be written as: 

6v g (r)=w^(r)-v g (r). (3) 

Similarly, one could derive an effective potential via 
an energy route, such that the full internal energy of a 
system a state point (iV, V, T), governed by the Hamilto- 
nian ([2]), is reproduced by the simpler two-body formula 

U(N,V,T) = X -p 2 J dridr 2 s(ri2)viKri2). (4) 

The difference with the bare-pair potential w^(r) can 
be written as: 

Svu{r) — (r) — vu(r). (5) 

It is not hard to show that both v g (r) and vjj (r) depend 
on the state point at which they are derived. Thus if 
they are used at a different state-point, one would expect 
transferability problems. What is perhaps more worrying 
is that one can also show that at a given state-point, v g (r) 
and vjj(r) cannot be the same. For example, To lowest 
order in p and it/ 3 -* , the ratio between the two corrections 

is [H, m, [SI: 

g$-i + ((.»)V). <e> 

Since v g (r) is unique, it is therefore impossible to repre- 
sent the properties of a system governed by the Hamil- 
tonian of Eq. ^ by a single pair potential. These rep- 
resentability problems are widespread in coarse-grained 
descriptions of soft-matter systems (29| . They are also 
important in other coarse-grained simulations. For ex- 
ample, one of us recently studied the structural route to 
derive radially symmetric potentials for water (36j. There 
v g (r) and vjj{r) were explicitly constructed and look very 
different, as anticipated by by Eq. [5]). At some relevant 
state-points, using v g (r) to calculate the virial pressure 
resulted in a compressibility factor Z that was almost 
two orders of magnitude larger than that of the origi- 
nal multi-site water model used to parameterise v g (r). 
Admittedly, it may not be surprising that an isotropic 



potential should perform so poorly when the underly- 
ing fluid has complex orientational correlations. On the 
other hand, one could derive a potential vp(r), designed 
to correctly predict the virial pressure of the underlying 
water model. But, due to the uniqueness of v g (r), this 
potential would no longer correctly reproduce the pair 
correlations. It is really a case of "choose your poison" . 

In practice, things may not always be so dire for many 
colloidal suspensions. First of all, one is not usually 
trying to reproduce all properties as accurately as one 
needs to for a molecular fluid like water. Secondly, for 
many hard-core colloids with short-ranged interactions, 
the state-dependence of the interactions is not expected 
to be too important. In this case, transferability and 
representability problems, which are often linked [29(, are 
not expected to be strong. 

For more complex soft-matter systems the story is 
more subtle, see e.g. a paper by Murtola et aZ.[37[ for 
an excellent discussion for the case of lipid bilayers. In 
general, careful coarse-graining often leads to potentials 
that are much softer than normal atomic or molecular 
fluids [l], 0, [H, Hi|. This can have important advan- 
tages when performing equilibrium calculations, because 
energy barriers are lowered, making, for example, Monte 
Carlo sampling more efficient. When used in dynami- 
cal simulations, however, the correct application of effec- 
tive potentials is more difficult to properly derive. For 
example, in the "polymers as soft-colloids" picture [HI, 
the inter-polymer interactions are so soft that polymers 
can easily pass through one another. This is clearly 
not physical, and so schemes that use effective poten- 
tials for polymer dynamics must re-introduce crossabil- 
ity constraints to prevent coarse-grained polymers from 
passing through one another [4Q] . Depletion interactions 
in a multi component solution also depend on the rela- 
tive rates of depletant diffusion coefficients and particle 
flow velocities |4ll . [42|. Thus care must be taken that 
the components integrated out can respond much more 
rapidly than the colloids move, so that an effective po- 
tential picture still approximately holds. 

In this context it is illuminating to consider the sta- 
tistical mechanical origins of dissipative particle d yan - 
mics (DPD), which also relies on soft interactions [43ll4^| . 
As discussed in paper I, this method includes two sep- 
arate innovations. The first is to use soft-potentials. 
As described above, these can indeed be very useful for 
speeding up sampling of equilibrium behaviour. For non- 
equilibrium behaviour however, the application of these 
potentials remains somewhat obscure. One hopes that 
the significant lowering of energy barriers leads nonethe- 
less to relative time-scales that are qualitatively correct. 
If this can indeed be shown to be the case, then the soft 
potentials are a very useful computational tool. In prac- 
tice however, this is often hard to show with any cer- 
tainty, which often limits the reliability or physical rele- 
vance of many DPD simulations. 

The second independent innovation in DPD arises 
from the use of a thermostat that conserves momentum. 
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Molecular dynamics simulations are most naturally per- 
formed in the microcanonical ensemble, but most investi- 
gators prefer to use other ensembles such as the canonical 
ensemble. To achieve this, they implement a thermostat. 
The problem is that many thermostats don't exactly con- 
serve momentum, and so the hydrodynamic interactions 
that would naturally be present in the microcanonical 
calculations are not rendered correctly. The DPD ther- 
mostat solves this problem. However, it can also be used 
with other inter-particle interactions, and doesn't depend 
on innovation 1 of DPD, the use of soft potentials. 

For colloids in solution however, using the DPD ther- 
mostat does not generate the correct hydrodynamics 
without a complete simulation of the underlying solvent. 
As pointed out before, this is highly inefficient, and so 
a different approach must be used, and that will be the 
subject of the next section. 



B. Coarse-graining dynamics of colloids 

To describe the dynamics of colloidal suspensions, we 
follow paper I and implement the SRD method first de- 
rived by Malevanets and Kapral [IH . The method works 
by exploiting the fact that Navier Stokes hydrodynam- 
ics arises from local momentum conservation. Thus one 
can employ greatly simplified dynamics that are compu- 
tationally very efficient to simulate, and still generate hy- 
drodynamic behaviour. There are many methods in this 
class, including direct simulation Monte Carlo (DSMC) 
method of Bird [HI, H(| and the Lattice Boltzmann (LB) 
technique where a linearized and pre-averaged Boltz- 
mann equation is discretised and solved on a lattice [47| . 
In particular Ladd and others have extended this method 

to model colloidal suspensions [am nam mm. 

SRD has the particular advantage that transport coef- 
ficients have been analytically calculated [Hj, HH, l5(| - 
greatly facilitating its use. It is important to remem- 
ber that for all these particle based methods, the par- 
ticles should not be viewed as some kind of composite 
supramolecular fluid units, but rather as coarse-grained 
Navier Stokes solvers (with noise in the case of SRD) [12] • 
The SRD fluid is modelled by N point particles of mass 
m, with positions rj and velocities Vj. The coarse grain- 
ing procedure consists of two steps, streaming and col- 
lision. During the streaming step, the positions of the 
fluid particles are updated via 



v t {t + 8t c ) = ri(t)+Vi(t)8t c 



(7) 



In the collision step, the particles are split up into cells 
with sides of length ao, and their velocities are rotated 
around an angle a with respect to the cell centre of mass 
velocity, 

Vi(i + St c ) = v c . m ,i(t) + Tli{a) [vi(t) - v c . m ,i(t)] (8) 

where v c .m,i = Sj'*( TOV i)/ J2 3 m ^ s the centre of mass 
velocity of the particles the cell to which i belongs, (a) 



is the cell rotational matrix and 5t c is the interval be- 
tween collisions. The purpose of this collision step is 
to transfer momentum between the fluid particles while 
conserving the energy and momentum of each cell. The 
algorithm is efficient because direct interactions between 
the particles are not taken into account. 

Spherical colloids of mass M can be embedded in an 
SRD solvent using a Molecular Dynamics technique, as 
first shown by Malevanets and Kapral [3| , and studied in 
detail in reference I. We employ the following scheme: For 
the colloid-colloid interaction we use a standard steeply 
repulsive potential of the form: 



<Pcc(r) 



1 (r > 2 1 /24 CTcc) 



while the interaction between the colloid and the solvent 
is described by a similar, but less steep, potential: 

Mr) JMW-OHM) (r 

\0 (r>2V«a cs ) 

where a cc and a cs are the colloid-colloid and colloid- 
solvent collision diameters. We propagate the ensuing 
equations of motion with a Velocity Verlet algorithm us- 
ing a molecular dynamic time step At 



Ri(t + At) = Ri(t) + V l {t)At - 



2M 



(9) 



vs + ^) = vm + m+ ^ + M) ^ do) 

where Ri and Vj are the position and velocity of the 
colloid, and Fi the total force exerted on the colloid. 
Coupling the colloids in this way leads to slip boundary 
conditions. Stick boundary conditions can also be im- 
plemented [57| . but for qualitative behaviour, we don't 
expect there to be important differences. In parallel the 
velocities and positions of the SRD particles are streamed 
in the external potential given by the colloids and the ex- 
ternal walls and updated with the SRD rotation-collision 
step every time-step St c . We chose a cc — 4.3ao- The 
choice of time-steps, as well as a number of other tech- 
nical issues are described in more detail in ref I, and for 
the case of two-dimensional discs, in[58|. 

One key step in the implementation and interpretation 
of these kinds of coarse-grained dynamical simulations 
lies in the realisation that hydrodynamic phenomena de- 
pend primarily on a set of dimensionless numbers, and 
that it is important to keep these in the right regimes in 
order to simulate the correct physical behaviour. Typi- 
cal numbers include the Reynolds number Re which mea- 
sures the relative strength of viscous and inertial forces, 
and the Mach number Ma, defined as the ratio of a typical 
colloidal velocity to the velocity of sound c s of the un- 
derlying fluid. In a real colloidal system, the Re and Ma 
numbers are very small, typically on the order of 1CP 6 or 
less. However, as long as one keeps them less than about 
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10 _1 , the physical regime doesn't really change much, a 
fact that can be exploited in simulations. For example 
if the real Ma number of a colloidal suspension were re- 
produced, it would mean that millions of sound waves 
would need to be resolved for modest motion of the col- 
loids. But to be in the right physical regime, a sound 
velocity well separated from other velocities is all that is 
needed. Obviously having a slower sound speed greatly 
increases the efficiency of the simulation. The example 
of the Ma number can be extended to a whole range of 
other hydrodynamic numbers, as discussed in more detail 
in reference 

A related set of arguments can be made for the time- 
scales of a typical colloidal system. A colloid of radius 
1/j.m diffuses over its radius in a time of order several sec- 
onds. At the same time, the shortest physically relevant 
time-scale for the colloidal dynamics is the Fokker- Planck 
time Tpp, defined in [ll| as the time over which the col- 
loidal force-force correlation function decays. For a typi- 
cal colloid, this may be as small as 10~ 13 s. Other impor- 
tant time-scales include the sonic time-scale t cs = R/c s , 
on the order of 10 _10 s for a typical colloid. The kine- 
matic time-scale r v — R 2 jv which measures the time for 
vorticity, which diffuses with the kinematic viscosity u, 
to diffuse over one colloidal radius is also important, and 
linked to the Re number. In a real colloidal system, these 
time-scales are all separated from each other by many or- 
ders of magnitude. However, to be in the correct physical 
regime, one doesn't necessarily need to resolve all those 
time-scales. All that is really needed is a clear time-scale 
separation. 

Our strategy is then as follows: First, the relevant 
time-scales are identified. These are then telescoped down 
to a hierarchy which is compacted to maximise simula- 
tion efficiency, but sufficiently separated to correctly re- 
solve the underlying physical behaviour. We illustrate 
this process in Figure [TJ and discuss it in more detail in 
ref I. Although the particular example studied here con- 
cerns a colloid in suspension, we argue that many other 
coarse-graining methods make implicit use of some as- 
pect of telescoping down. A careful time-scale analysis is 
very important for the correct physical interpretation of 
the simulation results. 

Having sketched our simulation technique, we now turn 
to an application where diffusion, convection, and hydro- 
dynamic interactions all play a role. 



III. TAYLOR DISPERSION FOR SMALL 
SOLUTES 

When a fluid flows between two parallel plates, stick 
boundary conditions mean that the local fluid velocity is 
zero at the plate surfaces. The resulting Poiseuille flow 
has a velocity distribution of the form [59( : 
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FIG. 1: Telescoping down: The hierarchy of time-scales for 
a colloid (here the example taken is for a colloid of radius 1 
fim in H2O) is compressed in the coarse-grained simulations 
to a more manageable separation. As long as the physically 
important times are clearly separated, the simulation should 
still generate the correct physical picture. Once the simu- 
lations are completed, they can be related in more detail to 
particular experiments by telescoping back out to the relevant 
experimental time-scales. 




(11) 



FIG. 2: Taylor dispersion between parallel plates. The initial 
solute gets stretched to a paraboloid shaped plug. Diffusion, 
indicated by the vertical arrows, evens out the concentration 
profile leading to a wider plug. Taylor dispersion describes 
the diffusive spreading of the plug. 



where u is the mean velocity, 2L is the plate separation, 
and y is a coordinate measuring the distance from the 
centre of the pipe. Thus the maximum velocity uq — |u 
is at the centre of the flow profile. 

The basic idea behind Taylor dispersion is that a solute 
particle will diffuse between the streamlines. Its average 
position moves forward with the average flow velocity, 
but because some streamlines move faster than the aver- 
age, and some slower, the diffusion between streamlines 
adds an extra dispersive component relative to its aver- 
age motion in the axial direction. 

To get a better physical understanding of Taylor diffu- 
sion, we follow an argument from [(icj . depicted in figEJ 
For time scales less than tjj ~ L 2 /D, convection initially 
stretches the solute into a parabola. The solute at the 
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front end of the parabola leads the solute at the edges by 
a distance Ugt after a time t. Diffusion across the channel 
then smears the parabola into a plug in a time-scale to 
it takes the solute to reach the edges. Thus the parabola 
is smeared into a plug of width WtD ~ uoto = uqL 2 / D. 
We can describe this as the the plug taking a random 
step of size ~ WtD at each 'step' to- If we now repeat 
the process for different stripes of the plug; each convec- 
tively stretched then diffusively smeared, then after N 
steps, the plug will have evolved as 



(W 2 ) ~ NW tD 



ulL 2 
D 



(12) 



Accordingly, in addition to molecular diffusivity D, the 
solute can be seen to diffuse along the channel, with re- 
spect to the average flow, with an effective dispersion 
coefficient 



D* 



<L 2 
D 



(13) 



Note that the effective dispersion coefficient is inversely 
proportional to the molecular diffusion coefficient, which 
may seem counter intuitive at first glance. However, so- 
lute molecules with higher molecular diffusivities spend 
more time diffusing across the pipe and less time sam- 
pling particular velocity streamlines thus reducing their 
effective spreading. 

A more sophisticated calculation based on flow in a 
three dimensional pipe of radius L yields a prefactor (6lj | 



D* = 



u 2 L 2 
48L> 



(14) 



In previous work [631 ] we have calculated the exact value 
for the prefactor in two dimensions. By integrating the 
equations of motion over the plate separation and imple- 
menting the appropriate boundary conditions we found 
that 



D* 



2u 2 L 2 
105D ' 



(15) 



In this paper we will carry out most of our simulations 
and analysis in two dimensions, simply because simula- 
tions are easier, and the fundamental physics we are after 
should not be very dependent on the exact dimension we 
work in. 

The calculations above assume that the molecular dif- 
fusion in the axial direction is negligible. When this is not 
the case, a good approximation for the Taylor dispersion 
coefficient was first given by Aris [64| to be: 



D 



TA 



D* + D, 



(16) 



and so Dta is often called the Taylor-Aris dispersion 
coefficient. 

Taylor dispersion has been extensively used to study 
the diffusion coefficients of small molecular solutes [65j 
where D is typically 10 _8 cm 2 / s so that D* is expected 
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FIG. 3: Finite size effect correction to the flow profile. Col- 
loids are no longer able to velocity gradients near the walls 
resulting in an average velocity higher than that of the sol- 
vent. 



to be large, making it easy to measure. Experimental 
setups examining atomic Taylor dispersion use tubes of 
order a few mm in diameter and lengths many meters. 
In such setups, D can be measured with an accuracy of 
±1% [24]. For colloidal systems, where diffusion coef- 
ficients are much smaller, it may be possible to use a 
shorter tube, because the Taylor-Aris dispersion coeffi- 
cient would in fact be much larger. 



IV. TAYLOR DISPERSION OF FINITE SIZED 
COLLOIDAL DISCS 

As illustrated in Fig. IIII1 when the colloid radius be- 
comes comparable to L, the colloids can no longer sample 
the entire Poiseuille flow profile. This leads to a number 
of corrections [25]. Firstly, because the colloids are ex- 
cluded from the slower flow-lines, their average velocity 
is faster than the underlying flow of the fluid. This can 
easily be calculated for Poiseuille flow 

L - R < y 2 \ , 3/2 2 1 o 

l -j 2 ) dy= r\? > + i x -i x 

(17) 

where \ — R/L measures the size of the colloid compared 
to the distance between the plates. For infinitely small 
particles with x — the average velocity is u c — u, but as 
X increases, so does the average velocity, until it reaches 
a maximum of u c = |u at \ = 1- This size exclusion 
phenomenon can used to analyse the relative diameters of 
colloidal particles (hydrodynamic chromatography [66]). 

For finite x, the Taylor dispersion coefficient also un- 
dergoes corrections [25j. If the new boundary conditions 
are taken into account, then our calculation of Eq. lfTT 
incurs further corrections 1631: 
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(18) 



for 2d flows. The radial dispersion of the colloids is 
now effectively reduced as they are unable to sample the 
higher velocity gradients near the edges, reducing the ra- 
dial dispersivity and leading to a lower value of the Taylor 
dispersion coefficient. 

Note, however, that the calculations above assumed 
that 1, the colloids perfectly follow the stream- lines, and 
are not themselves affected by their proximity to the 




FIG. 4: Simulation snapshots at different time-interals of the 
Taylor diffusion of colloids in a flow field. 



plates, and 2, that the average flow profile is not itself 
perturbed by the presence of the colloids. For three di- 
mensions, Brenner & Gaydos [67T | used a moment analysis 
to show that extra terms enter into Eq.(fl7|) that reduce 
the numerical coefficients so that the actual flow velocity 
is lower than a calculation which ignores the wall effects. 
What happens to the average flow in two dimensions, and 
how this exactly affects the Taylor dispersion beyond that 
shown in Eq. (|18[) . is not known. 

To study the effect of finite particle size we performed 
computer simulations of colloids and heavier SRD tracer 
particles using the methodology described in papers I 
and [58| . Fig. [W] shows some typical snapshots of col- 
loidal dispersion driven by Taylor dispersion. We take 
care to make sure that the Re number is always less than 
about 0.1, and that other hydrodynamic numbers are also 
in the correct regime. We simulate with up to 60 colloids 
and for box lengths of 238i?. Averages are measured after 
the transients are deemed to have decayed away. 

One effect that the colloids can have on the flow is that 
they increase the viscosity, which in turn affects the av- 
erage velocity. This is shown in the top panel of Fig. [S] 
By fitting the measured profile for different colloid pack- 
ing fractions <f> we find that the viscosity of the mixture 
scales roughly as rj^ ~ (1 + 1.95</j)?yo. In three dimensions 
and in the bulk, the pre-factor in front of (j> is 2.5 [ToL [Tll| . 

In the bottom panel of Fig. we show the average ve- 
locity of colloidal particles as a function of the radius to 
pipe-size ratio These show a decrease of their veloc- 
ity compared to the ideal calculation of Eq. (|T7)l due to 
the effect of hydrodynamic interactions with the walls. 
We also tested the average velocity of SRD tracer par- 
ticles ("ideal colloids"). These have a larger mass than 
the other fluid point particles, and in addition, interact 
with the walls with an interaction of radius R. Othewise 
they are point-particles and participate in the SRD colli- 
sion rules. We expect that these particles have a smaller 
effect on the fluid flow profile than the colloids do. As 
shown in Fig. these particles also don't follow the ideal 
correction term, and so presumably also experience hy- 
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FIG. 5: Top : Effect of the colloid concentration on the sol- 
vent flow profile. The straight line represents the theoretical 
flow profile for a pure SRD solvent. The simulation points 
represent the flow profile for a homogeneous distribution of 
colloids at a packing fraction of <j> = O.f while the dotted line 
shows the fit with the increased viscosity. Bottom : The av- 
erage flow velocities of colloids (circles) and "ideal" colloids 
(squares) are compared to the original simple prediction of 
Eq. [17] that ignores the effects of the walls. As expected, for 
increasing \ — R/L, the walls on average slow down the par- 
ticles compared to what one would expect if they just followed 
the unperturbed Poiseuille flow lines. 



drodynamic interactions with the wall that slow them 
down. 

Within a simple Langevin picture, a colloid loses mem- 
ory its velocity on a time-scale given by te = M/£, where 
the friction coefficient £ is defined by D = ksT/^. Al- 
though for colloids this Langevin picture is in fact heavily 
misleading, see e.g. the discussion in the appendix of rcf 
lfl3|. for our purposes here it should suffice as a useful 
reference time-scale. 

If the fluctuation-dissipation theorem is invoked, then 
an effective Taylor dispersion friction £* can be defined, 
such that D* — ^r-. In a simple Langevin picture, fluc- 
tuations will die out on the time scale 



M 



MD* 2 MD u 2 L 2 



k B T 105 k B T D 2 105 
where we have defined the pipe Peclet number 
uL 

Pe = — = tu/T DL 
which can be viewed as the ratio of the time tj 



t^Pe 2 (19) 
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FIG. 6: Top panel: Integral of the velocity fluctuations. Bot- 
tom panel: Effective mean square displacement with respect 
to the average flow. Both graphs are for "ideal colloids" made 
up of heavier SRD particles. 




FIG. 7: Effective mean square displacement of colloids under- 
going Taylor dispersion. Simulations were carried out in pipes 
of sizes x = 0.358,0.179,0.134,0.0719,0.0448 respectively. 
The solid straight line represents the approximation for Tay- 
lor dispersion from Eq. (|15|l . and the dotted lines denote the 
finite size correction predicted in Eq. (|18[1 . Simulation param- 
eters were chosen to ensure that the self and Taylor dispersion 
are always clearly separated, i.e. that D*/D = tg* /t% > 5. 
All runs were carried out for colloid Reynolds number of less 
than 0.1. 



it takes the fluid to move in the axial direction over a 
distance L, to the time tdl it takes a particle to diffuse 
that distance. 

To measure Taylor dispersion as accurately as possi- 
ble, we simulated 'ideal' colloids of mass M = 62.3m 
under flow and measured the dispersion of 25 of these in 
a pipe of width such that x = 0.134. A separate set of 
simulations was done in the same pipe, but with "ideal" 
SRD colloids with the same wall-colloid interaction as 
the solvent. For these ideal particles and the relatively 
wide pipe, we don't expect strong wall effects, so that the 
simple Eq. (|18p should provide a good approximation to 
the Taylor dispersion. The results of these simulations 
are shown in FigJBJ The plain curves denote results for 
the heavier SRD fluid particles without an extra wall- 
interaction, while the dashed curves depict results for 
"ideal colloids" , i.e. the heavier SRD particles with an 
extra interaction with the wall. In the top section of 
the graph, we have plotted the integral of the velocity 
auto-correlation function for both sets of particles. The 
horizontal lines are a guide to the eye and intersect the 
y-abscissas at the theoretical values for dispersion pre- 
dicted in Eq. (Tlg|) for \ — and x = 0.134 respectively. 
In both instances, we see that the simulated integrals 
plateau near the lines suggesting good agreement with 
Taylor theory, including the corrections of Eq. (fT5|) . The 
lower part of Fig[5] shows the mean squared displacement 
of the particles w.r.t. the average of the centre of mass of 



a large number of colloids. This is expected to scale as 
X 2 ~ 2D*t. Note that even for this relatively small value 
of Xi the Taylor dispersion coefficient is down by about 
|, so that the effect of a finite size is much stronger on 
the Taylor dispersion than it is on the relative velocity of 
the particles, which only shows a very modest correction 
for that value of %. 

Fig [7| depicts the mean square displacement w.r.t. the 
centre of mass flow of normal colloidal discs. The mea- 
sured Taylor dispersion coefficient D* appears to lie be- 
tween the predictions of Eq. (fT5|) and (|18p . 

We also investigated the behaviour of the velocity au- 
tocorrelation function(VACF) C(t) = (6v(t)dv(Q)} for 
colloids undergoing Taylor dispersion in a pipe of width 
X = 0.13437. These are shown in Fig[5J and compared to 
the VACF for particles for the case of no flow. For short 
times both curves are well described by a simple Enskog 
kinetic equation (see ref [HI]). Physically this is to be 
expected, as these kinetic processes occur on the sonic 
time-scale which is much faster than the flow, and so it 
is gratifying to see that our simulations correctly sepa- 
rate out this time-scale from the others. At long times 
the two VACFs differ. As shown in the bottom panel, the 
particle without flow shows the expected t -1 long time 
tail, first observed by Alder and Wainwright 68] . This 
tail is caused by hydrodynamic correlations that prop- 
agate through the fluid. For a finite sized pipe, we ex- 
pect the tail to feel the influence of the walls at longer 



9 




FIG. 8: Decay of the velocity autocorrelation function for 
particles undergoing normal (circles) and Taylor (squares) dis- 
persion. The solid line in the top panel depicts the Enskog 
exponential decay. The dotted line in the bottom panel shows 
t -1 decay and serves as a guide to the eye. 

times [H, HI], and some evidence of this may be visi- 
ble. Nevertheless, this tail clearly differs from that of 
the VACF for Taylor dispersion, which shows a much 
more rapid non-algebraic decay. Instead it decays on a 
time scale comparable to . The reason for this is most 
likely that the hydrodynamic correlations are washed out 
by the effect of the flow. 



than the flow when the ratio \ = R/L increases, but 
that hydrodynamic wall effects also slow down the col- 
loids when the pipes become very narrow. Finite x nas 
an even bigger effect on the Taylor dispersion coefficient, 
decreasing its value from that expected for \ — 0- We 
also measured the VACF under flow conditions, and find 
that this decays more rapidly than the case of no flow, 
most likely because the Taylor dispersion process breaks 
up the long-time hydrodynamic correlations. 

Finally, given that we used a "telescoping down" of 
the hierarchy of time-scales to derive a tractable simu- 
lation scheme for the dynamics of our colloidal system, 
how would one "telescope up" to make experimentally 
relevant predictions? The first step is to cast properties 
of interest into dimensionless form as much as possible. 
From that we see, for example, that the relative decrease 
of Taylor dispersion depends only on the ratio %, and 
not on the absolute value of either the pipe size or the 
colloid radius R. Thus many different experimental col- 
loid and pipe sizes should be describable by the same 
set of simulations. On the other hand, if we are inter- 
ested in the long-time tails, then a more subtle analysis 
is needed. It appears that for colloids undergoing flow, 
the VACF decays on a time of the order of t^*, which, 
through Eq. (fT9|) , can be related to the Langevin time 
and the Pe number. Again, one can rather straightfor- 
wardly estimate for colloids if one knows their diffusion 
coefficient. However, the relative values of t% and the dif- 
fusion time td that a colloid needs to diffuse over its 
own radius may be very different in the simulation than 
in an experiment. That means that the two time-scales 
in the simulation need to each be scaled differently when 
comparing to experiment. So in the end, to interpret a 
coarse-grained simulation, the best advice is to try and 
understand as best as possible the underlying physics. 



V. SUMMARY 

In summary then, we have adapted a hybrid MD-SRD 
technique to study colloidal discs in two dimensions un- 
dergoing flow in a narrow pipe. By carefully consider- 
ing the time and length-scales involved in the problem, 
we were able to derive a coarse-graining method that is 
tractable, but still resolves the main physical properties 
under investigation. We measured corrections to the av- 
erage flow profile, and find that the colloids move faster 
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